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We perform a comparative study of the free energies and the density distributions in hard sphere 
crystals using Monte Carlo simulations and density functional theory (employing Fundamental Mea- 
sure functionals). Using a recently introduced technique (Schilling and Schmid, J. Chem. Phys 131, 
£S) ' 231102 (2009)) we obtain crystal free energies to a high precision. The free energies from Fun- 

damental Measure theory are in good agreement with the simulation results and demonstrate the 
applicability of these functionals to the treatment of other problems involving crystallization. The 
agreement between FMT and simulations on the level of the free energies is also reflected in the 
density distributions around single lattice sites. Overall, the peak widths and anisotropy signs for 
C**) ■ different lattice directions agree, however, it is found that Fundamental Measure theory gives slightly 

narrower peaks with more anisotropy than seen in the simulations. Among the three types of Fun- 
damental Measure functionals studied, only the White Bear II functional (Hansen-Goos and Roth, 
J. Phys.: Condens. Matter 18, 8413 (2006)) exhibits sensible results for the equilibrium vacancy 
^ ' concentration and a physical behavior of the chemical potential in crystals constrained by a fixed 

vacancy concentration. 

g PACS numbers: 82.70Dd,61.50Ah,71.15Mb 

■ I. INTRODUCTION 

O ' 

The phase behavior of hard spheres is one of the most intensely studied subjects within the realm of classical 
statistical mechanics. The existence of a fluid-solid transition has been predicted already more than fifty years ago by 
early computer simulation methods [l|, 0] ■ Advances in colloidal engineering have led to the experimental realization 
of almost hard sphere-like systems, confirming the occurence of crystallization in such systems in the 1980's Q. The 
variety of hard sphere-like colloidal systems include polymeric spheres (index-matched solvent, sterically stabilized) 
[|| and also thermotropic colloids @ using particles with diameters of the order of a few hundred nm. This allows 

■ the use of scattering techniques with visible light and/or the use of real-space microscopy to resolve single-particle 
| positions. Using these systems and techniques, numerous features of the statics and dynamics of the crystallization 
. process and the competing glass transition have been studied in detail (see e.g. Refs. |6h10(). 

The progress in real-space imaging opens the perspective that the static density distribution in crystals and the 
dynamics of the nucleation process can be studied with unprecedented resolution. The primary information obtained 
in these experiments, the trajectories of single particles, is very much the same as the information obtained in a 
computer simulation. Thus the further analysis of this primary information brings together these two fields. Currently, 
' e.g. the processes of homogeneous and heterogeneous nucleation in colloidal hard sphere systems are under scrutiny 
[111 1 1 21 ] , and the unambiguous resolution of the underlying mechanisms of these processes appears to be possible using 
simulation/real-space experiment on the one side and the established reciprocal-space (scattering) experiments on 
the other. 

From the theory side, classical density functional theory (DFT) is a good candidate to study crystallization phe- 
nomena on a microscopic level. The concepts of equilibrium DFT have been developed over the past forty years (for 
an early review see Ref. [HI)- In this context, hard spheres appear to be one of the few classical fluids for which 
quantitatively predictive density functionals can be constructed thanks to powerful geometric arguments, leading 
to the so-called Fundamental Measure Theory (for recent reviews see Refs. [2 EH)- In contrast to the maturity 
of equilibrium DFT, dynamic DFT is a still developing field which has has been started only about ten years ago 
[F6l - |2(j |. Centerpiece of dynamic DFT is the time evolution of the inhomogeneous one-particle density. The most 
intensely studied variant of the theory is actually an approximation to Brownian dynamics, thus it appears to be 
well-suited for the study of colloidal systems. However, due to the complexity of the FMT functionals, any dynamic 
DFT studies of a hard sphere system with inhomogencities in two or three dimensions have not been undertaken. 
In fact, there are only a few equilibrium studies of inhomogeneous problems in two and three dimensions [2lTl24j ]. 
unrelated to the crystallization problem. The study of hard sphere crystals within FMT has been restricted so far to 
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sensible paramctrizations of the density distribution in a cr ysta l, nevertheless this approach has elucidated the key 
features of a reliable DFT for the crystallization transition J25||. (A more detailed review of the problem of crystal 
phases within density functional theory is given below.) 

In order to make progress in the direction of applying dynamic DFT (with the FMT functionals that work very well 
for hard spheres) to the currently studied nucleation problems (TTl Il2j. we will study first the more modest problem 
of the static density distribution in hard sphere crystals in this paper. This will be done by a full, three-dimensional 
minimization of the FMT functionals and contrasted to the results of our Monte-Carlo simulations. (Surprisingly, 
the density distribution in hard sphere crystals has been likewise studied very little using simulations.) Such a study 
is an absolute prerequisite for the more difficult dynamic problems involving crystal-fluid interfaces to be tackled in 
the future. We will show that the full minimization discriminates between different FMT functionals which are very 
similar in the description of the fluid phase. We will shed some new light on the problem of an equilibrium vacancy 
concentration within DFT. We will demonstrate that the FMT results for the free energy per particle for the crystal 
phase are in very good agreement with the corresponding simulation result which has been produced by a recently 
introduced method. 

The paper is structured as follows. In Sec. |n] we briefly review the density functional approach to crystallization 
in the hard sphere system. Sec. IIIII discusses a few points relevant for the FMT crystal description in more depth. 
These address the constrained minimization in the unit cell (with particle number fixed) , the relation of the respective 
constrained chemical potential to Widom's trick in a system with fixed vacancy concentration, and the numerical 
procedure of the FMT functional minimization. In Sec. HVI we briefly describe our Monte Carlo method to obtain free 
energies and density distributions. Sec. |V] compiles our results on free energies, equilibrium vacancy concentrations 
and density distributions and in Sec. I VII we present our conclusions. 

II. HARD SPHERE CRYSTALS IN DENSITY FUNCTIONAL THEORY 

In density functional theory, the crystal is viewed as a self-sustained inhomogeneous fluid, i.e. an inhomogeneous 
density profile p C r(r) minimizes the grand potential functional 

n[>(r)] = T\p(T)\ - J d 3 rp{v) {fx - V oxt (r)) , (1) 

with the external potential V ext being zero. Here, p is the chemical potential and F[p] is the free energy functional 
which is conventionally split into an ideal and an excess part: 

F\p] = F A [p\+^[p] (2) 

with the exact form of the ideal part given by 

F {d [p] = J d 3 r/ id (r) = j d 3 rp(r)(hx[p(r)A 3 ]-l) . (3) 

Here, A is the de-Broglic wavelength. It was realized very early (in 1979) that a simple Taylor-expanded version of 
the excess free energy 

= J d 3 r/ ex (r) « ~J d 3 rj dVc$(r - r'; p rof )Ap(r)Ap(r') (4) 

allows for minimizing solutions p C r(r) (26l |. Here, p ro f is a reference density around the liquid coexistence density 
and c rc{ is the direct correlation function in the bulk liquid at this reference density (which is related to the bulk 

(2) 

structure factor by S(k) = l/p Ie f — l/c£efW)' I n this early work, the minimization to obtain p cr was a constrained 
one: expanding the density as 

Pcr(r) = Po+Y^ Pi ex P( iK j ' r ) ( 5 ) 

3 

(Kj is the set of the reciprocal lattice vectors) , the minimization was only performed with respect to the moments pj 
which belong to the first or to the first and fourth shell of reciprocal lattice vectors. 1 In this approximation, one sees 



1 Reciprocal lattice vectors belong to the same shell if they transform into each other under the point group transformations from the 
considered crystal symmetry. The first shell contains all reciprocal lattice vectors with the lowest magnitude, etc. As an example, for 
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that the crystal free energy in Eq. (|4| "probes" the Fourier transform c^](k) only at one or two values of k. For fee 
these values are k\ « 10.9/a and k± 20.8/a where a is the side length of the cubic unit cell, very near the first two 

maxima of c r J(k) resp. S(k). At first sight, it may appear surprising that an expansion of the free energy like Eq. (@|, 
valid at small density variations is sufficient to sustain the rapidly varying density profile in a crystal. However, 
the isotropic correlations between two particles in Fourier space are described by the structure factor (and hence by 

c). ei ). Since the shells of reciprocal lattice vectors for an fee lattice are also distributed fairly isotropically, the possible 
description of a solid with a density near the reference density appears to be less unexpected. Subsequent work has 
revealed that the expansion in reciprocal space ([S]) is converging slowly. Furthermore, there are serious quantitative 
problems in this approach if it comes to the description of the lattice density peak width (much too narrow), crystals 
at higher density (unstable) and the vacancy density (around 10 percent at coexistence which is a factor of about 100 
too large) [13, El- 

A more general approach to inhomogeneous hard sphere fluids in general and to the description of crystals in 
particular consists in the ansatz 



F™[p] = J d J rp(r)*|/(p(r)) . (6) 

Here, \P is a suitable function of a weighted density 

p(r) = J d 3 r'p(r')w(r — r'; p) = p * w (r) (7) 

which employs a weight function which in turn may depend on the weighted density itself. (For the Taylor expanded 
functional ((¥]), p — const. + c). c{ * p.) The functions \& and w can be determined through the equation of state and 
the bulk direct correlation functions which are assumed to be known. Here, the self-consistent solution for w may 
be rather involved in particular realizations. Examples for this class of functionals include the Tarazona functionals 
Mark I [H] and Mark II [30|, the weighted-density approximation (WDA) [H[ and the modified WDA [13]. Crystal 
structures in these approaches have been usually obtained by minimizing the ansatz 

3 

N ( — ) 2 exp (— a(r — r,) 2 ) (8) 



lattice sites i 



with respect to the Gaussian peak width a and the normalization N. If rt vac denotes the relative concentration of 
vacancies then N = 1 — 7i vac . With such a Gaussian ansatz, the reciprocal lattice modes of the density (see Eq. ((3J)) 
are given by pj = N exp(— K 2 /(4a)). Using the most sophisticated versions of these weighted-density approaches, 
one can achieve a rather good agreement with simulations for the liquid-solid coexisting densities and a physically 
sensible behavior also for denser crystals. This is understandable since in comparison with the simple Taylor expanded 
functional (jj) the weighted-density form © includes contributions from higher-order direct correlation functions 
and through the self-consistent determination of w it is guaranteed that at higher densities the changed isotropic 
correlations in a (possibly metastable) reference liquid are taken into account. Still, the lattice density peaks come 
out too narrow compared to simulations, the crystal free energy per particle is too small by about 5% and a quasifree 
minimization in modified WDA (with the restriction n vac = 0) revealed qualitatively wrong peak asymmetries in 
the different lattice directions as well as an unphysically large interstitial density [331 ] . A sensible, small vacancy 
concentration n vaCi o which minimizes the free e nerg y can only be obtained by incorporating an appropriate additional 
constraint term into the free energy functional [34| . 

However, the WDA approach which is built on the isotropic fluid correlations cannot be expected to treat coordina- 
tion effects in crystals correctly on a fundamental level. These include the description of the metastable hard sphere 
bec crystal [35[ and the crystal-fluid interface [36|. Here, the development of fundamental measure theory (FMT) 
marks an important breakthrough [37| . FMT postulates an excess free energy with a local free energy density in a 
set of weighted densities n a : 

F cx [p] = J d 3 r<S>(n a (r)) . (9) 



fee, the reciprocal lattice is bee and the first shell contains 8 reciprocal lattice vectors 27r/a(±f, ±f , ±f) where a is the side length of 
the cubic unit cell. The fourth shell (used in Ref. [3) contains 24 reciprocal lattice vectors, given by 27r/a(±3, ±1,±1) plus two cyclic 
permutations of the Cartesian components. 
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FIG. 1: Three types of cavities which can hold only one particle. 



The weighted densities are again constructed as convolutions of the density with weight functions, n a (r) = p*w a (r). 
In contrast to the WDA, the weight functions reflect the geometric properties of the individual hard spheres (and not 
the properties of an interacting pair). For one species, the weight functions include four scalar functions w° . . .u> 3 , 
two vector functions w^w 1 and a tensor function w 1 defined as 

2 2 

w 3 = 9{R - |r|) , iv 2 = 6{R- |r|) , w 1 = ^— , w° - W 



4ttR ' AirR 2 
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w 2 = j^5(R- |r|) , w : = 



r| v 1 " ' 4ttR 



w io=^ S (R-\r\). (10) 



Here, R is the hard sphere radius. Using these weight functions, corresponding scalar weighted densities no ■■■ri3, 
vector weighted densities ri! , n 2 and one tensor weighted density n t are defined. In constructing the free energy density 
arguments concerning the correlations in the bulk fluid and arguments for strongly inhomogeneous systems are 
used. For the bulk, $ is required to reproduce exactly the second and third virial coefficent of the direct correlation 
function. Furthermore, imposition of the Carnahan-Starling equation of state (38l . l39j and/or consistency with a 
scaled particle argument [13, H(| leads to a closed form for $. The arguments using strongly inhomogeneous systems 
are known in the literature under the label "dimensional crossover" [4l|: Through a suitable external potential, the 
hard sphere fluid can be confined to lower dimensions and the density functionals for these lower-dimensional systems 
should emerge from the correct density functional in 3d. Of particular relevance are the crossover to Id where the 



exact density functional is known [42[ and to Od where a hard s phe re is confined to a point and the free energy is 



a simple function of the mean occupation number at this point [4JJ . The Od confinement can be realized through 
differently shaped cavities (overlapping spheres of radius R) which can hold only one particle (see Fig.[l|. Respecting 
the Od limit for different cavities is of particular relevance for the crystal description since this means that the mutual 
exclusion of hard spheres in various coordinations is correctly described. In Refs. (25l . [43j a solution is given which 
respects the Od limit for cavities (a) and (b) of Fig. Q] and approximates the Od limit for cavity (c). 

The arguments presented in the above paragraph lead to the following form of the excess free energy density: 

$({n[p(r)]}) = -n m(l - n 3 ) + tpi(n 3 ) " 



1 - n 3 

3 (—712 n 2 • n 2 + n2,int,ijn2,j + n 2 n tjij n t ji — nt,ijnt,jknt,ki) 
16tt(1 - n 3 ) 2 



¥> 2 (n 3 ) — ! — ~ ^ ■ ( n ) 



Here, (fii(n 3 ) and ^2(^3) are functions of the local packing density n 3 (r). With the choice 

<pi = l, <P2 = l (12) 



we obtain the Tarazona tensor functional [25( which is built upon the original Rosenfeld functional 37f . The latter 
gives the fluid equation of state and pair structure of the Percus-Yevick approximation. Upon setting 



<pi = 1 (13) 
-2n 3 + 3n 2 3 - 2(1 - n 3 ) 2 ln(l - n 3 ) 

4 



^ = 1 3nT 
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we obtain the tensor version of the White Bear functional |38|, consistent with the quasi-exact Carnahan-Starling 
equation of state. Finally, with 

2n 3 -Tt§ + 2(l-n3)ln(l-n 3 ) . . 

tpx — 1 H (14) 

3n 3 

_^ 2n 3 - 3nj + 2n| + 2(1 - rc 3 ) 2 ln(l - n 3 ) 

the tensor version of the recently introduced White Bear II functiona l 1401 is recovered. This functional is most 
consistent with restrictions imposed by morphological thermodynamics |44j . 



III. HARD SPHERE CRYSTALS IN FMT 
A. Minimization and p consistency 

The minimization of the grand potential functional 



Q[p] = T id [p] + F^[p] - J d 3 r(p - V^(r))p(r) , (15) 

leads to 

/T 1 ln(p oq (r)A 3 ) = ~p cx [ Pcq (r)} +fx- V^(r) . (16) 
The functional p[p(r)] is given by 

^[P(r)} = -^j^ (17) 
Sp{r) 



"-'E/*'^ W ° {T '- Th (18) 



In principle, for a force-free system (V ext = 0), the specification of a suitable chemical potential p should lead upon 
minimization to a periodic crystal profile p cq = p C r(j) with the bulk density po(//)- The side length a of the cubic 
unit cell and consequently the vacancy concentration n vac should adjust itself to comply with Eq. (|16p . Here, n vac is 
connected to the occupation of the unit cell of the fee lattice by 

d 3 rp(r) = 4(l-n vac ) . (19) 

coll 

In practice, such a procedure is not feasible. Rather, for a given bulk density po, also n vac (and thus a) is prescribed 
and a constrained free energy functional for the unit cell 

^'lceil= / d 3 rf id [p}+ f d\r[p]-p' ( d s r(p(v)- Po ) (20) 

./cell ./cell ./cell 

is minimized where p! = p! (po,n vac ) plays the role of a Lagrange multiplier to ensure (|19[) . In the work reviewed 
previously // was not determined explicitly, and only a few studies bothered to vary also n vac (which should be 
close to zero) such that the free energy per particle is indeed minimized. However, there is a useful consistency 
condition between p.' and p(po)- Let / cr (po) denote the free energy density for the fully minimized crystal with 
vacancy concentration n vaC: o- Then 

M = ^ = M'(Po,"vac)|„ vac= „ vac0 (21) 

This can be shown as follows. Let p eq (po, n vac ; r) be the minimizing density profile for a crystal with fixed bulk density 
Pa and vacancy concentration n vac . Using the expansion in reciprocal lattice vectors ([5]), p = po + Ylj Pj cxp(iKj • r), 
it is seen that the constrained minimization of Eq. (|20|) yields 



M/9o,n vac ) = r / dr 

a d (P0,«vac) ./cell d P0 

1 



a 3 (P0,^vac) icell 



d 6 r(Hp eq A d )+p ex [ Peq ]) . (22) 
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Here, the last line follows since J d 3 r df ex /dpo = J d 3 r (SJ 70 * /5p) (dp/dpo) and dp/dpo = 1. On the other hand, the 
chemical potential from the crystal equation of state becomes: 



dfcr df c[ dn 
dpo dn vac dp 



= Pq- 



n vac — n vaCi o 

d(F CI /N) On, 
dn m , 



d 3 r 5(/ id [Pcr]+/ CX [Pcr]) 



dpo 



a 3 (Po,"vac,o) Jcoii dPo 

+ A*'(P0, ™vac,o) = fJ-'(PO, ™vac,o) 



(23) 



n vac — n vaCj o 



Thus we see that p(fo) = P-'ipo, n v&c,o)y since the crystal free energy particle per particle (F CT /N) is minimal at 

^vac ^vac.O- 



B. Basic considerations on single defects 



As we have seen in the previous considerations, the appearance of defects enters the equilibrium density profile in a 
crystal through the average occupation of a lattice site. The dominating type of defect in the equilibrium hard sphere 
crystal arc monovacancies whose properties have been studied before in simulations explicitly [45l - l47j . In order to 
derive a general formula for the constrained chemical potential p'(po-, ^vac,o) it is useful to discuss the thermodynamics 
of a crystal containing vacancies more in detail. 

Here we follow Ref. [46| in the subsequent reasoning. We introduce a system with M lattice sites which contains 
n monovacancies at given positions and index thermodynamic quantities with these numbers, such that e.g. Vm,i 
denotes the volume of a lattice with M sites, 1 fixed vacancy and therefore M — 1 particles. Furthermore it is 
convenient to define by -/ vac the change in free energy due to the creation of a single vacancy at a specific lattice 
point while keeping the volume and the number of lattice sites constant: 

- /vac - F M+hl (M,V M +ifi,T) - F M+lfi (M + l,V M +i,o,T) 

= -m(p A 3 )-/ v c a x c . (24) 

As usual, the free energy F(N 7 V, T) is a function of particle number N , volume and temperature. In the second line 
we have separated — / vac into the ideal gas contribution and the excess part —/vac- Assuming no interaction between 
pairs of monovacancies, the total free energy Fm,u is 

Fm.u = F M ,o - "/vac = Mf - n/ vac (25) 

where /o is the free energy per particle (or per lattice site) in a defect-free crystal (i.e. precisely the value of F/N 
determined in our simulations). In order to calculate the equilibrium concentration of vacancies, it is more convenient 
to switch to the Gibbs free energy Gm,u{M — n,p,T) in a system of M — n particles at constant pressure p and 
temperature T. We define g vac as the change in G due to the creation of a single vacancy at a specific lattice point: 

.9vac = Gm+i,i(M,p,T) — Gm,o(M,p,T) 

= F M +i,i{M,V M +i,i,T) - F M>0 (M,VMfi,T)+p(V M+ i,i - V u ,o) ■ (26) 

Using Eq. (J2H) and furthermore / = F M +i,o(M + 1, V M +i,o,T) - F M fi(M, V Mfi ,T) and p, a = f Q +pV M ,o/M we find 

5vac = MO — /vac ■ (27) 

The total Gibbs free energy G^ n includes the entropic contribution due to the distribution of n vacancies over M 
lattice sites (n <C M): 

G\f n « G M -„,o + ngvac + nk B T (in ^ - l) (28) 

Minimizing with respect to n yields the equilibrium concentration of monovacancies n vac : 

n vac .o = jj = cxp(-/?g vac ) = exp(-/?(/i - / vac )) . (29) 

Let us now define an excess chemical potential p.' wi (n vac ) for a constrained crystal at a fixed vacancy concentration 
n vac through the free energy of particle insertion (Widom's trick). Its excess part can be estimated by the probability 
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Pace (Vws) °f inserting a particle into the Wigner-Seitz cell (with volume Vws) around the vacancy position and the 
probability n vac of picking the vacancy lattice site among all lattice sites. Thus: 

Mwi( n vac) ~ ln(p A 3 ) - k B T\n(P acc (Vws)) - fc B Tlnn vac 

= /vac - fc B Tlnn vac (30) 

The second line follows since — k B T hi P acc (Vws) is precisely the excess free energy cost /° ac of removal of one vacancy 
[461 ]. We see immediately that in equilibrium, n vac = n vac .o, we have Pwi( n vac,o) ~ Mo which demonstrates the con- 
sistency between the thermodynamic and insertion route in equilibrium. (Note that the correction to the equilibrium 



chemical potential is only linear in n vac .o [461] .) However, for the constrained system the insertion route predicts that 
//(po, ™vac) diverges upon n vac 0. 

The system with the constraint of fixed n va c corresponds to the free energy functional in Eq. (|20j) and thus we 
may identify Mwi( n vac) = A*'(Po> "-vac)- Therefore the logarithmic increase of the chemical potential with n vac — > 
is a stringent test for fully minimized density functional models. However, we want to point out that physically the 
divergence of // with vanishing vacancy density is not entirely correct as outlined in the following. Even in a perfect 
lattice it is possible to insert another interstitial particle. Similarly to — / vac one can define the change in free energy 
/in due to the creation of a single interstitial at a specific lattice point while keeping the volume and the number of 
lattice sites constant: 



/in = F MA (M + 1, V MA ,T) - F M , (M, V M ,a,T) . (31) 

The second index for F and V refers to the number of interstitial particles in the system. Therefore it follows that 
for vanishing vacancy concentration the constrained chemical potential is given by 

/AP0,™vac < ™vac,o) = /in + 0{ ^vac ) • 

(32) 

Simulation results for the free energies / vac and /j n in hard sphere crystals near coexistence give approximately the 
magnitudes 8 k B T and 34 k B T, respectively O, HtJ. Since f- m S> f vac , it is clear that the constrained chemical 
potential should exhibit the logarithmic divergence upon n vac — > down to very small vacancy concentrations. (At 
coexistence, p'(po, "vac.o) = Mo ~ 16 k B T. For smaller n vac , \j! should rise up to approximately f m ps 34 k B T and 
then level off.) 

C. Previous results in FMT 

A fully three-dimensional minimization of FMT aiming at the crystal profile has not been carried out before. In 
Tarazona's ground-breaking work [25j the density profile was parametrized as 

Pcr(r) = (l-«vac) (^expHr-r,) 2 ) (l + AV/^r-r,)) , (33) 

lattice sites i 

U(r = (x, Vl z)) = .t 4 + + z 4 - j|r 4 . (34) 

Here, fa is the leading term for the unit cell anisotropy in cubic lattices. The free energy per particle was minimized 
with respect to n vac , a and K4 using the Rosenfeld tensor functional (TiTj) and (fl"2"]) . The anisotropy turned out to 
be unimportant for the values of F CT /N (modifying it by less than 10 -3 k B T). Both F CT /N and a were shown to be 
in good agreement with the old simulation data of Ref. (48j . No clear free energy minimum was found for a nonzero 
n vac , indicating that n vac .o < 10~ 8 . 

Concerning the issue of the equilibrium vacancy concentration n vaCj o in FMT, there are two more, partially contra- 
dictory statements in the literature. In Ref. [4l| it was argued that the correct Od limit of a density functional (for 
particles strictly localized to their lattice sites) should always lead to a finite, but small n vac .o. The Od excess free 
energy is given by (3F^ = 77+ (1 — 77) ln(l— 77) with a corresponding excess chemical potential /3/i§5 = — ln(l — 77). Since 
the packing fraction at each lattice site is corresponds to 1 — ?i vac , one finds /?Mod = — l nn vac- The equilibrium vacancy 
concentration follows upon identification of /igd( n vac,o) with the crystal chemical potential /io as n va c.o = ex p( — P^ad)- 
According to this argument one would expect an equilibrium vacancy concentration n vac .o ~ 10 -8 at coexistence. 
We observe that the divergence of /3po5 is precisely of the type derived before for the constrained chemical potential 
//(po, rivac) (see Eq. ([BUI) ). However, the plain identification /j^j = /j' is incorrect due to the neglect of the free energy 
of vacancy formation. This explains the four orders of magnitude difference in rt vac .o when compared with simulations 
f4|||47j]. Another approach was taken in Ref. [49| to calculate n vaCj o- There, the Rosenfeld functional (among others) 
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was minimized in a perturbative approach assuming isotropic density distributions around lattice sites and an expan- 
sion around the close-packing limit. A free energy minimum was found for values of 7i. va c,o consistent with simulation 
results. However, we will demonstrate below that this finding is not consistent with our full minimizations. 

The success of the density parametrization using isotropic Gaussians and zero vacancy concentration inspired the 
works of Ref. [50j to investigate non-fcc crystals and of Refs. [HI H2 to treat binary systems and the crystal-fluid 
interface within FMT. In the latter work, the interface density profile was parametrized in an intuitive way, however, 
in this way one cannot ensure that crystal and fluid are in chemical equilibrium (see Sec. IVl below). 

D. Numerical solution of the FMT Euler— Lagrange equation 

In actual calculations, we determine the constrained crystal profile p C q(po, nvac] r) by a full minimization in 
three-dimensional real space. For such a three-dimensional problem, the density profile p and 11 weighted den- 
sities (two scalar densities ri2, n^, three vector densities (ri2)i for i = {x,y, z} and six tensor densities (jit)ij for 
{ij} = {xx, yy, zz, xy, xz, yz}) need to be discrctized on a three-dimensional grid covering the cubic unit cell. Usually 
we chose grids with dimensions 64 3 (for lower densities around the coexistence density PcoexC 3 ~ 1.04) up to 256 3 
(for higher densities). Here, a = 2R is the hard sphere diameter. The necessary convolutions were computed using 
Fast Fourier Transforms. With prescribed po and n vac , the constrained functional (|20[) is minimized through Picard 
iteration (with mixing) of the Euler-Lagrange equation (|16[) . A new profile pi+\ is determined from an old profile pi 
and an appropriate p! i through 

p l+1 = ap' i+1 + (1 - a) pi , (35) 
Pi+i = CX P ( -Z^^ylA] + ) • (36) 

Here, u\ is determined such that J" u d 3 r pi+i = 4(1 — n vac ). The mixing parameter a is of the order of 0.01. The 
iteration was stopped when the relative deviation between p' t and p! (po,n vac ) from Eq. (|22[) was below 5 ■ 10~ 6 . The 
iteration procedure was stabilized by two means: (i) enforcing the physical requirement n 3 (r) < 1 — n vac at each 
iteration step since the singularity at = 1 (sec the functional in Eq. (1111) ) is avoided in that manner, (ii) enforcing 
the point symmetry of the fee crystal in the density profile in each iteration step. In each iteration, this point symmetry 
is slightly violated by numerical inaccuracies. Without correction and using a ~ 0.01, this symmetry violation quickly 
grows and eventually leads to numerical singularities. Only for a < 10~ 5 , convergence was achieved without explicit 
enforcement of the point symmetry at the price of an increase in computation time by a factor 100-1000. For a given 
Po, n vac is varied and the location of the minimum is checked using the consistency condition ([21) . However, for 
^vac < 10~ 5 it proved to be hard to arrive at a convergent solution. 

IV. MONTE CARLO SIMULATIONS 
A. Computation of absolute free energies 

The computation of absolute free energies poses a problem to Monte Carlo simulation, because it requires the 
evaluation of the partition function. For most systems that have an infinite and continuous state space, the partition 
function cannot be computed directly. However, one can compute free energy differences and derivatives by MC 
simulation. Hence, if there is a suitable reference system, the free energy of which is known analytically, free energies 
can be extracted from MC simulation. Here we use a technique that was recently introduced by Schilling and Schmid 
[55| . The technique extends the well-established thermodynamic integration with respect to the harmonic crystal 
(Einstein crystal) [HH to disordered reference states and tether potentials that are not harmonic. We compare our 
results to DFT and to simulation results obtained by Vega and Noya using the Einstein Molecule (EM) technique (a 
variant of the Einstein crystal that avoids having to correct for the center of mass motion of the system [HI). 

We used systems of size N = 1728 and potential wells of radius r cuto ff = 0.75 a. The path of the thermodynamic 
integration for each density was subdivided into simulations of 35 different values of the coupling strength e between 
and 80, each of which consisted of 10 4 equilibration sweeps and 10 6 sweeps of averaging (where one sweep consisted of 
N attempted particle moves). There are two sources of error: The first one results from the prodecure of integration 
and can be estimated to (j3AF)i nt /N — 0.001. The second is statistical. These errors were computated by using 
the Jackknife algorithm with 1000 subsets on each part of every integration. With this the statistical errors equally 
amount to (pAF) sta . t /N = 0.001. 
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TABLE I: Comparison of free energies (a) calculated using the algorithm of [55| and (b) the results obtained with the EM 
method from [5J]. For the DFT results, the White Bear II functional and n vac = 1CF 4 was used. 'Gauss' refers to minimization 
using the Gaussian approximation (Eq. (H)). The free energies according to the Speedy equation of state have been determined 
using Eq. (55) ■ In order to obtain numbers, A = a has been used. 



Poo" 3 


R /77/\/"( a ) 


pr j iv 

( \j ori/1 

[1\ — zU4o J 


pr j iV 

yl\ — > OO J 


R ZTV^^ / AT 

(Gauss) 


R ZTU^^ / AT 

(full min.) 


1.00 


4.530(2) 






4.541 


4.539 


1.04086 


4.960(2) 


4.955(1) 


4.9590(2) 


4.979 


4.977 


1.049 


5.048(2) 






5.069 


5.067 


1.08 


5.397(2) 






5.424 


5.422 


1.09975 


5.631(2) 


5.627(1) 


5.631(1) 


5.660 


5.658 


1.11 


5.756(2) 






5.787 


5.785 


1.14 


6.142(2) 






6.174 


6.172 


1.15000 


6.277(2) 


6.269(1) 


6.273(2) 


6.310 


6.308 



PF Spccdy /N 

4.532 
4.961 
5.049 
5.398 
5.631 
5.756 
6.140 
6.275 



B. Density Distribution 



A second set of simulations has been carried out to sample the density distribution of the hard sphere crystal unit 
cell with high accuracy. A perfect fee crystal has been set up in a cubic box with fixed side lengths, periodic boundary 
conditions and a particle number of N = 4 • n 3 Particles, which corresponds to n unit cells along one side of the box. 

The simulations have been run on a standard octocore CPU. The results for the density distributions presented below 
are based on simulations of a system consisting of N n= i3 = 8788 particles. Different system sizes of N n —g = 2916, 
N n= n = 5324 or iV n= is = 13500 particles have been used to study the extrapolation of the average Gaussian width 
a (see Eq. (j8|)) for N oo. Snapshots of the system configuration have been taken every 30 sweeps. After each 
sweep, appropriate global shifts of particle coordinates were applied to keep the center of mass fixed. The simulated 
bulk densities vary from p^a 3 = 1.04 to p^a 3 = 1.30, with the corresponding density distributions averaged over 
(1 . . .2) • 10 11 snapshots (depending on exact acceptance rate). This corresponds to approximately 18 days of CPU 
time for a single density distribution. Error estimates for po& 3 = 1.04 and pq<j 3 = 1.20 have been obtained using the 
Jackknife algorithm, with the unit cell histograms divided into 1000 subsets. All simulations have been started with 
100,000 sweeps of equilibration. In order to obtain the density distributions in a single unit cell, the system snapshots 
of all unit cells were mapped onto one unit cell providing us with a 3D histogram with a resolution of 80 bins per unit 
cell length. 



V. RESULTS 



A. Free Energies and coexistence densities 



In order to connect to the previous simulation work in Ref. [54j ]. we studied hard sphere systems with particle 
densities po<J 3 = 1.04086, 1.09975 and 1.15. Additionally, we also considered more density points and the values for 
the free energy per particle F/N for all our calculations are reported in Tab. Q] All simulation results were obtained 
with n vac = 0, whereas in the DFT results (using the White Bear II (WBII) functional) n vac = 10~ 4 was chosen 
which did not affect the values for F/N to the accuracy shown. Our simulation results and the results from Ref. [54| 
are consistent with each other on the level of 0.05 %. The DFT results are systematically larger than the simulation 
results, the discrepancy here is also not larger than 0.5%. 

The last column in Tab. U gives the corresponding free energy results as derived from the popular equation of state 
proposed by Speedy (56j . This equation of state is given in the form 

7) c 2 

a ^ , (37) 

(38) 



/3pSpcody 


/3Psing 


P 


P 




3 


P 


Pep 1 
n 1 
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TABLE II: Coexisting fluid (pa) and crystal (p C r) densities (the corresponding packing fractions are given in brackets), as well 
as the chemical potential /x coox and the pressure p C oox at coexistence for the three investigated DFT models. Here, RF is the 
tensor modified Rosenfeld functional with the free energy density determined by Eqs. (JTTJ and (|12[) , WB is the tensor modified 
White Bear functional (Eqs. |TT} and fLi)) ) , and WBII is the tensor modified White Bear II functional (Eqs. (JTTJ and CLU). 
The MC results are taken from Ref. [5q[. In order to obtain numbers, A = a has been used. 





PflO" 3 (7/fl) 


pcrO" 3 (??cr) 


/3/-£cocx 




RF 


0.892 (0.467) 


0.984 (0.515) 


14.42 


9.92 


WB 


0.934 (0.489) 


1.022 (0.535) 


15.75 


11.28 


WBII 


0.945 (0.495) 


1.040 (0.544) 


16.40 


11.89 


MC 


0.940 (0.492) 


1.041 (0.545) 




11.576 



where p cp cr 3 = V2 is the close-packing density, and c\ = 0.5914, ci = 0.7079 and C3 = 0.6022 are fitting constants 
determined recently from a fit to a substantial set of pressure data [57| . The pressure pspeedy exhibits the divergent be- 
havior (psing) for p — > p cp as predict by free-volume considerations. In order to obtain the free energy, thermodynamic 
integration can be applied after the divergent piece has been subtracted and integrated separately: 

eedy / •> f P /3(PSpeedy ~ Psing) 

'Pep 



{p) = , ^ r*°M - 3 ln[( Pcp - p)o*] + C . (39) 



Here, C is an integration constant which has been quoted in the literature [59[ from a fairly old simulation [60( 
as C — 2.843 ± 0.040. Fitting C by using Eq. (|39)) to our MC free energy data we find the improved estimate 
C = 2.8247 ±0.0006. 

For DFT, we used the Gaussian approximation to the density profiles to calculate the thermodynamic properties 
of liquid-solid coexistence via the Maxwell construction. The result is given in Tab. |TT] for the three investigated 
fundamental measure models and compared to very recent simulation results [58l|. F or the tensor modified Rosenfeld 
(RF) and White Bear (WB) functionals we recover the results quoted in Refs. |25l. |38| . For the RF functional, the 
coexistence densities are substantially smaller than the corresponding densities for WB and WBII. This is entirely due 
to the insufficient accuracy of the Percus-Yevick equation of state on the fluid side which underlies the RF functional. 
The WB and WBII functionals reduce to the Carnahan-Starling equation of state for homogeneous densities and thus 
their thermodynamic description of liquid-solid coexistence is satisfactory. 



B. Vacancy concentration and constrained chemical potential 

In Sec. 1111 Bl we have derived an expression for the constrained chemical potential p'(po,n vac ) (see Eq. (|30[) L 
Furthermore we recall that p'(pQ,n vac ) is precisely the Lagrange multiplier in the constrained minimization of the 
unit cell free energy, see Eq. (|20|) . We have examined its dependence on n vac and the bulk density p for the three 
functionals with the surprising result that only in the case of the White Bear II functional p'(po, n vac ) shows a weakly 
divergent behavior as n vac — > 0. The divergence appears to be weaker than — lnn vac , however. Furthermore only for 
the White Bear II functional the consistency condition (|2"Tj) is fulfilled for a (small and) finite equilibrium vacancy 
concentration n V ac,o- There is no minimum for the free energy per particle F/N upon variation of n vac for the cases of 
the Rosenfeld and the White Bear functional (neither in the Gaussian approximation nor for full minimization). This is 
consistent with Tarazona's finding of no minimum for n vac > 10~ 8 using the Gaussian approximation in the Rosenfeld 
functional (25[. As an exemplary result, we show p'(po,n vac ) for poc 3 = 1-04 (coexistence) for the three functionals, 
see Fig. [2] (a). The large discrepancies between results for the three functionals is somewhat surprising, given the 
fact that F/N varies only very little (in the Gaussian approximation we have (3F/N = 4.929 [RF], 4.912 [WB], 4.970 
[WBII] at this density). Our results for the explicit minimization also question the reliability of the approach taken 
in Ref. [49| to calculate n vac ^. It appears that the equilibrium vacancy concentrations and corresponding free energy 
minima are artefacts of the approximations used therein (isotropic density distributions around lattice sites and an 
expansion around the close-packing limit). 

In Fig. [5] (b) we show the variation of the equilibrium vacancy concentration with the bulk density for the WBII 
functional. There is reasonable agreement between the Gaussian approximation and the full minimization. However, 
the predicted n vac .o is consistently smaller (up to one order of magnitude) than available simulation results. Never- 
theless one should keep in mind that the DFT results do not follow from an explicit computation of the free energy 
of a vacancy / vac (see Eq. (f2"4"])) as the simulations do. It would be interesting in the future to calculate / vac through 
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FIG. 2: (a): The constrained chemical potential p'(po, ?i vac ) as obtained by full minimization of the three DFT models at 
the bulk density pocr 3 = 1.04. The dashed line shows the value for the chemical potential following from the thermodynamic 
definition of p = df CI /dpo where / cr is the free energy density. It is equal for the three DFT models to the given accuracy. If and 
only if p'(po, Hvac) = p, the free energy per particle is minimal and thus thermodynamic consistency holds, (b): Equilibrium 
vacancy concentration vs. bulk density as obtained for the WBII functional (full line-Gaussian approximation, filled diamonds- 
full minimization) and compared to available simulation results (open diamonds-Ref. [47|], filled squares-Ref. [HI). 



an explicit minimization of DFT around a fixed vacancy. Note that an initial attempt in that direction has been 
undertaken in Ref. [6l[ using the MWDA. 

An important implication arises from the fact that the consistency condition (|21[) , p'(po, w V ac,o) = df CI /dpo, can be 
fulfilled only for the WBII functional. It means that a free DFT minimization of the fluid-crystal interface which is 
consistent with the coexistence data from the Maxwell construction (see Tab. [TTJ) will not be possible with the WB 
and the RF functional . This follows since the fluid chemical potential at coexistence does not match the crystal 
chemical potential obtained by full mimimization. 



C. Density distributions 



The density distribution in the hard sphere crystal consists of nearly isolated density peaks around the lattice sites, 

Pcr(r)= ]T p(r-r 4 ) (40) 

lattice sites i 

with no appreciable overlap in the tails of p(r). In first approximation, p(r) is a Gaussian with a width parameter a, 

3 

p(r) ps p G {r) = (-) ' cxp(-ar 2 ) . (41) 



We will analyze the deviations from the Gaussian form in terms of an average radial deviation /agM an d an 
anisotropic deviation / a niso(r): 

p(r) ps PG {r) f AG (r) / aniso (r) . (42) 

The average radial deviation will be parametrized as 

f AG (r)=ex V [b 2 ar 2 +b i (ar 2 ) 2 + b 6 (ar 2 ) 3 ] , (43) 

where 62,^4,^6 ^ 1 are expected to be small. For the analysis of the directional anisotropy we apply a polynomial 
expansion in the form: 

/ aniso (r) = 1 + K, a 2 (x* + y 4 + z 4 - \A + K 6 a 3 (x e + if + z 6 - V) . (44) 
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FIG. 3: (a): Logarithm of the Gaussian width parameter a vs. bulk density po: DFT-WBII in Gaussian approximation (full 
line), DFT-WBII in full minimization (circles), extrapolation to the thermodynamic limit in MC (+ symbols) and results from 
Ref. [4^1 (squares), (b): radial probability r 2 p(r) in [100] direction for the bulk density pocr 3 = 1.04. Comparison between DFT 
and our simulations. 



This corresponds to the leading two terms in the cubic cell asymmetry (consistent with the point symmetry of the 
fee lattice). 2 



1. Gaussian width parameter 



For the DFT results, the width parameter a is the only minimization parameter in the Gaussian approximation 
once the normalization is fixed. From the results of the full minimization we determined a by a global fit with the 
Gaussian form (|4T|) to the lattice peak density distribution. The same was done using the MC data, additionally the 
value Q!oo in the thermodynamic limit was determined by the extrapolation from the values at finite box length L 
through the relation cxn = A/N 1 ^ 3 + cvoo [48}. For the densities po<j 3 = 1.05 and 1.13 we checked and confirmed this 
scaling for the four values N = 2916, 5324, 8788 and 13500. For the other bulk density values, we used N = 5324 and 
13500 to determine aoo- 

In Fig. [3] (a) we compare the DFT results for a with from our MC simulations and a corresponding width 
parameter extracted from the work of Young and Alder [48| . There the mean square deviation was determined which 
we converted to the Gaussian width parameter by assuming the Gaussian form for p(r): a = 3/(2(r 2 )). There is 
excellent agreement between the two simulations and also fair agreement between DFT and simulations. The Gaussian 
peaks in DFT are narrower than the simulated peaks which is similar to (M)WDA results although in (M)WDA the 
quantitative deviation is already considerable (compare e.g. with Tab. I in Ref. [62]). The radial probability, 
proportional to r 2 p(r), along the [100] direction is shown in Fig. [3] (b). As a remark, earlier MC data for the radial 
probability were erroneously scaled in the graphical presentations of Ref. (33| . 

In order to quantify th e sp read of the density distribution around a solid peak it is convenient to define the 
Lindemann parameter [63l |64| as the dimensionless root mean-square displacement: 



C = —J f d 3 rr 2 p(y) . (45) 

r nn y JwSC 



2 The density distribution around a lattice site can be expanded as p(r) = po( r ) + Pi( r )%i + 2~Zij Pij (rfxiXj + . . . , where Xi = Xi/r 
and the expansion coefficients Pij...( r ) are isotropic functions. From symmetry we have p(r) = p(— r) and py...(r) = pp(i)p(j)...( r ) 
where P is a permutation of the Cartesian indices. This implies that all expansion coefficients with an odd number of indices are 
zero and that p\\ = pi2 = P33, giving only an isotropic correction to second order (which can be absorbed into poM)- The lowest 
nontrivial expansion coefficients arc pim and P1122 which are not independent of each other since 1 = (x\ + x\ + X3) 2 . In our 
fits, we have chosen the radial dependence pmi(r) = K^a 2 pQ (r) /ag i r ) r ^ ar[ d also demanded that the angular integral of the 
anisotropy corrections over the unit sphere vanishes. This leads to an isotropic offset, such that the isotropic piece for our case becomes 
Po(r) = PG (r) /agM (1 - (3/5)K 4 « 2 r 4 ). 
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FIG. 4: Lindemann parameter £ vs. bulk density po for Monte Carlo simulation (MC) and DFT, also in comparison with 
Ref. [33[ | (Ohnesorge et. al.). Data for (M)WDA from Ref. [3^1 are taken for full minimization. 

Here the spatial integration is over a Wigner-Seitz cell (WSC) centered around a lattice position at the origin 

1 /3 

and r nn ~ o (v^/po) denotes the distance between two nearest neighbours in the crystal lattice. Data for the 
Lindemann parameter £ versus bulk density po are shown in Fig. 2) Clearly, £ is about 0.13 at melting and decreases 
with increasing density. The Monte Carlo data published earlier in Ref. [33j agree with those from our simulations. 
All density functionals considered here (WDA, MWDA, WBII) yield Lindemann parameters which are only slightly 
lower than the simulation data. However, the density profiles from (M)WDA and WBII differ, the agreement in the 
value of £ is due to an unphysical high interstitial density in the (M)WDA profiles |33j |. 



2. Deviations from the Gaussian form and anisotropy 



In Figs. [5] and IH1 we show in an exemplary way the density distributions in the principal lattice directions [100], 
[110] and [111] for the densities p^a 3 = 1.04 (near coexistence) and poa 3 = 1.20, respectively. The simulation data are 
always very close to the Gaussian form with the coefficient 64 of the leading deviation from the Gaussian form being 
small, [64! < 0.01 (see panel (a) in Figs. [5] and IB]) ■ Interestingly, b 4 changes sign at around p^a 3 = 1.10, indicating that 
below that density the distribution is wider than a Gaussian (larger curtosis) and above that density the distribution 
is narrower than a Gaussian (smaller curtosis). In DFT- WBII the density distribution has a smaller curtosis than 
a Gaussian with 64 « —0.03 for the range of densities 1.04 to 1.20 (see panel (c) in Figs. [5] and [5]). Turning to the 
asymmetries we note that our ansatz (Eq. (0D)) describes the data for small distances r very well and starts to deviate 
only when the overall density has dropped by a factor 10 4 compared to the center of the peak. This is illustrated 
in panels (b) and (d) in Figs. [5] and [6] where we compare the fit to the anisotropic part / an iso to the quotient of the 
density profile with the purely radial fit, p{v) / {pc{r) /agW ) (see Eq. (|42j) ). The qualitative behavior of the density 
distribution in the principal lattice directions is the same for MC and DFT- WBII, only the magnitude of the leading 
anisotropy coefficient A"4 is larger in DFT-WBII by about a factor 1.7. The agreement in sign and order of magnitude 
in A4 with simulations distinguishes fundamental measure theory from the (M)WDA approach where an opposite 
sign is obtained (33[. (Intuitively, the density distribution in [110] direction should be narrower than in [100] since in 
[110] direction the next neighbor is closer.) 

In Fig. [7] we show the value of A4 for a range of bulk densities from the fits to both MC and DFT-WBII results. 
The scatter in the data is a result of the uncertainty in the fits, but one can clearly observe a trend to lower K4 
for higher density. This would be consistent with the observation in Ref. [48| that towards close-packing the density 
distribution becomes Gaussian (A4 = 0). We note furthermore that we could not extract any meaningful results for 
the next-to-leading anisotropy coefficient K§ whose modulus appears to be smaller than IAT4I but the error estimate 
is always of about the same magnitude. Finally, we remark that our results for K4 are in quantitative agreement to 
earlier computer simulation data published in Ref. [33| . 





FIG. 5: Lattice site density distributions along the lattice directions [100], [110] and [111] for the bulk density poo -3 = 1.04. 
Panels (a) and (b) show MC results (N = 8788), panels (c) and (d) results from DFT-WBII. Panels (a) and (c) show p vs. r 2 in 
logarithmic scale, thus illustrating the deviation from a Gaussian form (straight line). The full line here is a fit to the Gaussian 
form pa (Eq. (gT|) with the parameter a = 77.5 (MC) and a = 84.4 (DFT-WBII). The dashed line is a fit to the non-Gaussian 
form pg/ag (see Eq. (g3j) with the parameters b 2 = -0.011, fe 4 = 0.0021, b 6 = -0.0002 (MC) and b 2 = 0.090, fe 4 = -0.029, 
&6 = 0.0009 (DFT-WBII). Panels (b) and (d) show the density along the three lattice directions divided by pg/ag- The 
lines show the corresponding anisotropics along the three lattice directions resulting from a fit to /aniso (see Eq. (04])) with the 
parameter K 4 = 0.022 (MC) and K 4 = 0.039 (DFT-WBII). 



VI. DISCUSSION AND CONCLUSION 



In this work we have performed a comparative study of the free energies and the density distributions in hard sphere 
crystals using Monte Carlo simulations and density functional theory (employing Fundamental Measure functionals). 
Using a recently introduced simulation technique, we could obtain crystal free energies to a high precision (see Tab. [XJ) 
which are consistent with the most recent parametrizations of empirical equations of state and allowed us to determine 
the crystal free energy in the close-packing limit with a higher accuracy than before (see Eq. ([39]) ). The free energies 
from Fundamental Measure theory are also in good agreement with the simulation results and demonstrate the 
applicability of these functionals to the treatment of other problems involving crystallization. The agreement between 
FMT and simulations on the level of the free energies is also reflected in the density distributions around single lattice 
sites (see Figs. [5] and [6]). Overall, the peak widths and anisotropy signs for different lattice directions agree, it is found 
that FMT gives slightly narrower peaks with more anisotropy than seen in the simulations. 

The deviations we observe between simulation and FMT point to possibilities of further improvement in the FMT 
functionals. Tarazona's construction of the tensor part of these functionals is an approximate representation of the 
three-cavity overlap situation (see Fig. [T]) which leads to a complicated expression. It would be interesting to study 
the close-packing limit of this expression in a systematic manner. 

Additionally we studied theoretically for the constrained minimization in the unit cell (with particle number fixed) 
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FIG. 6: Lattice site density distributions along the lattice directions [100], [110] and [111] for the bulk density pou 3 = 1.20. 
Panels (a) and (b) show MC results (N — 8788), panels (c) and (d) results from DFT-WBII. Panels (a) and (c) show p vs. 
r' 2 in logarithmic scale, thus illustrating the deviation from a Gaussian form (straight line). The full line here is a fit to the 
Gaussian form p G (Eq. (glj) with the parameter a = 343.7 (MC) and a = 399.0 (DFT-WBII). The dashed line is a fit to the 
non-Gaussian form p G /ag (see Eq. with the parameters b 2 = 0.014, b 4 = -0.0054, 6 6 = -0.00002 (MC) and b 2 = 0.075, 

6 4 = -0.026 , 6 6 = -0.0002 (DFT-WBII). Panels (b) and (d) show the density along the three lattice directions divided by 
Pg/ag- The lines show the corresponding anisotropies along the three lattice directions resulting from a fit to / an iso (see 
Eq. (gU)) with the parameter K 4 = 0.014 (MC) and K 4 = 0.025 (DFT-WBII). 

the relation of the respective constrained chemical potential /j,' to Widom's trick in a system with fixed vacancy 
concentration n vac . The latter analysis gives a simple relation, fx' = const. — lnn vac (see Eq. (|30p ). which poses a 
consistency condition on the corresponding FMT results for It turns out that from the three studied variants of 
FMT, only the White Bear II functional shows the qualitatively correct behavior whereas the Rosenfeld and the White 
Bear functional give qualitatively incorrect results (see Fig. [2]). This implies that for further studies such as the free 
minimization of the crystal-fluid interface or nuclcation processes only the White Bear II functional is a promising 
candidate. 
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